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ABSTRACT 

This paper presents two families of phase-space distribution functions that generate 
scale-free spheroidal mass densities in scale- free spherical potentials. The assumption 
of a spherical potential has the advantage that all integrals of motion are known explic- 
itly. The 'case F distribution functions are anisotropic generalizations of the flattened 
f{E, Lz) model, which they include as a special case. The 'case IF distribution func- 
tions generate flattened constant-anisotropy models. Free parameters control the radial 
power-law slopes of the mass density and potential, the flattening of the mass distribu- 
tion, and the velocity dispersion anisotropy. The models can describe the outer parts 
of galaxies and the density cusp structure near a central black hole, but also provide 
general insight into the dynamical properties of flattened systems. Because of their sim- 
plicity they provide a useful complementary approach to the construction of flattened 
self-consistent three-integral models for elliptical galaxies. 

The dependence of the intrinsic and projected properties on the model parameters 
and the inclination is described. The case I models have a larger ratio of rms tangen- 
tial to radial motion in the equatorial plane than on the symmetry axis, the more so 
for smaller axial ratios. The case II models have a constant ratio of rms tangential 
to radial motion throughout the system, as characterized by Binney's parameter p. 
The maximum possible ratio fp/cTp of the mean projected line-of-sight velocity and 
velocity dispersion on the projected major axis always decreases with increasing radial 
anisotropy. The observed ratio of the rms projected line-of-sight velocities on the pro- 
jected major and minor axes of elliptical galaxies is best fit by the case II models with 
f3 ^ 0. These models also predict non-Gaussian velocity profile shapes consistent with 
existing observations. 

The distribution functions are used to model the galaxies NGC 2434 (El) and 
NGC 3706 (E4), for which stellar kinematical measurements out to two effective radii 
indicate the presence of dark halos (CaroUo et al.). The velocity profile shapes of both 
galaxies can be well fit by radially anisotropic case II models with a spherical logarithmic 
potential. This contrasts with the f{E,Lz) models studied previously, which require 
flattened dark halos to fit the data. 

Key words: galaxies: elliptical and lenticular, cD - galaxies: individual: NGC 2434 
and NGC 3706 - galaxies: kinematics and dynamics - galaxies: structure - line: profiles. 



1 INTRODUCTION 

Elliptical galaxies are dynamically complex systems (e.g., de Zeeuw & Franx 1991). Many unsolved problems regarding their 
structure still exist. Studies of the presence and properties of dark halos and massive central black holes have been hampered 
by lack of information about the stellar velocity dispersion anisotropy. However, the body of observational data from which 
such knowledge can be derived (at least in principle) is growing steadily. In particular, deviations of the shapes of the stellar 
line-of-sight velocity distributions, or 'velocity profiles' (VPs), from Gaussians can now be measured reliably (e.g., van der 
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Marel et al. 1994a; Bender, Saglia & Gerhard 1994; CaroUo et al. 1995). Detailed dynamical models based on phase-space 
distribution functions (DFs) are needed to interpret such data. 

Very few, if any, galaxies are spherical. For a proper interpretation of the high-quality data that can now be obtained, 
flattened axisyrnrnctric models appear a minimum requirement. Axisymmetric potentials admit two classical integrals of 
motion, the energy per unit mass E, and the angular momentum component per unit mass parallel to the symmetry axis, L^. 
In addition, there is usually a non-classical third integral, /a (Binney & Trcrnainc 1987). The DF / is generally a function 
f{E, Lzjis). A subset of all DFs is formed by those that depend only on the two classical integrals, i.e., f{E, Lz). Construction 
of such models has only recently become practical, and as a result they are now being studied in groat detail (e.g.. Hunter & 
Qian 1993; Evans 1993, 1994; Dehnen & Gerhard 1994; Evans & de Zeeuw 1994; Qian et al. 1995; Dehnen 1995; Kuijkcn 1995; 
Magorrian 1995). These models properly take flattening into account, and are very useful for their relative simplicity. However, 
they are special in that the radial and vertical velocity dispersions are everywhere equal (i.e., an = az)- This is generally not 
the case in real galaxies (Binney, Davies & lUingworth 1990; van der Marel 1991). For most practical applications (e.g., to 
demonstrate that no model without a black hole or dark halo can fit a particular set of observations) one needs to construct 
more general models with three-integral DFs. This is complicated, because I3 is generally not known analytically, and may 
not even exist for all orbits. Dejonghe & de Zeeuw (1988) considered axisymmetric models with a potential of Stackel form, in 
which a global third integral is known analytically. Dehnen & Gerhard (1993) considered axisymmetric models with a flattened 
isochrone potential in which the third integral can be approximated using Hamiltonian perturbation theory. Models of this 
kind are currently being applied to fit the kinematics of real galaxies (e.g., Dejonghe et al. 1995; Matthias & Gerhard 1995). An 
alternative to these semi-analytical methods is a completely numerical approach in which individual orbits are superimposed 
(as in Schwarzschild 1979) using a linear programming, maximum entropy, non-negative least squares, or related technique. 
Several groups are now working on this. All these methods require substantial analytic or numerical efl^ort. 

A case of 'intermediate complexity' which has not been much studied in the literature is that of flattened mass distributions 
in spherical potentials. The assumption of a spherical potential has the advantage that the squared angular momentum per 
unit mass, L^, is an explicitly known third integral, so that / = f{E,L^,Lz). White (1985) already determined simple DF 
components in the spherical logarithmic potential (see also de Zeeuw, Evans & Schwarzschild 1995), while Kochanek (1994) 
solved the Jeans equations for realistic flattened density distributions in this potential. Mathieu, Dejonghe & Hui (1995) 
constructed triaxial mass models in a spherical potential for the galaxy Cen A. In the present paper a detailed study is 
presented of DFs for scale-free axisymmetric mass densities embedded in scale-free spherical potentials. The general form 
of the DFs of such models is derived, and two particular families of DFs are discussed in detail. For these families most 
physically and observationally interesting quantities can be determined analytically. The models therefore allow a detailed 
study of the dependence of the observable kinematical quantities on the various model parameters, such as the power-law 
slopes of the mass density and potential, the axial ratio of the density distribution, the inclination angle of the symmetry axis 
with respect to the line of sight, the intrinsic velocity (lispc>rsiou aiiisotropy, and the amount of mean streaming. The models 
provide significant insight into the dynamical structure of flattened galaxies, and provide a useful complementary approach 
to the construction of fully self-consistent models. 

As an application, the issue of the evidence for massive dark halos around elliptical galaxies is considered. Both tangential 
anisotropy and the presence of a dark halo can cause the observed velocity dispersion to remain roughly constant out to well 
beyond the effective radius. Hence, a flat observed dispersion profile does not prove the existence of a dark halo. Carollo et 
al. (1995) presented observations of the major axis kinematics and VP shapes for four elliptical galaxies, and constructed 
f{E, Lz) models to interpret their data. Here we restrict the discussion to two of the four galaxies, NGC 2434 and 3706, for 
which Carollo et al. concluded from the observed VP shapes that dark halos must be present. They showed that the dark 
halos must be flattened, if the observed VPs are to be fit with an / = f{E, Lz) DF. We use our models to determine whether 
the dark halos of these galaxies must indeed be flattened, or whether the VPs can be fit as well with a spherical dark halo 
and a DF of the form f{E, L^, Lz). 

Section 2 discusses the DFs of the models, and the calculation of the intrinsic and projected velocity moments and VP 
shape parameters. This section is technical, and readers mostly interested in applications of the models might wish to skip 
to Section 3 after Section 2.1. Section 3 gives a general discussion of the various properties of the models. The application to 
the Carollo et al. data is described in Section 4. Section 5 summarizes the main results. Appendix A presents an algorithm 
for the reconstruction of a VP from its moments. 



2 THE MODELS 

2.1 Potential and mass density 

Throughout this paper {r,9,(j)) denote the usual spherical coordinates, and {R,(j),z) the usual cylindrical coordinates, with 2 
along the symmetry axis of the mass density. The relative potential ^ and the relative energy per unit mass £ = — ^v^ are 
defined as in Binney & Tremaine (1987). The potential decreases outwards. Its value at infinity can either be finite, in which 
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case one can set = without loss of generality, or it can be *I'oo = —oo. The quantity £ is the binding energy per unit 
mass of a star. Only stars with £ > ^too are bound to the system. 
Scale-free spherical potentials ^' are considered, of the form 

M/M-V^ X / -ln(r/b), (5 = 0; 

nr) = Vo X I ^^^^^^ (1) 

where & is a reference length. The free parameter S, with < S < 1, determines the radial slope of the potential. The scale 
velocity Vo is equal to the circular velocity at the reference length. For 5 = the potential is the logarithmic potential, and 
the circular velocity is independent of radius. For 5 = 1 the potential is Keplerian. If the Kepler potential is generated by a 
total mass M, then Vq = GM/b, with G the universal constant of gravitation. 
Mass densities are considered that are power laws on oblate spheroids: 

p(7?,.) = po(f + ^)--''^ (2) 

where po is a reference mass density, 7 > is a constant that determines the radial fall off, and g < 1 is the constant axial 
ratio of the similar concentric isodensity surfaces of the mass distribution. The eccentricity is e = -^/l — q^. The limit g = 1, 
or e = 0, describes the spherical power law p{r) = po Mass distributions of the form (^) always produce systems with 

infinite total mass: for < 7 < 3, the total mass diverges at large radii, for 7 = 3, the total mass diverges at both small and 
large radii, while for 7 > 3, the total mass diverges at small radii. Nonetheless, the models meaningfully describe the properties 
of realistic finite-mass systems, but only at those radii where the mass density can be approximated by equation (^. 

Dimensionless quantities are used throughout the remainder of this paper: f = r/b; R = R/b; z = z/b; v = u/\/2Vb; 
L = L/V^bVo; £ = £/V(f; * = */V^; p = p/po; and / = f / po{2Vo)~^^'^ ■ Henceforth, aU tildes are omitted. The potential 
and mass density of the models are thus: 

T / \ I ~ ln»", 5 = 0: 

"^'■^^i^r'^ 0<S<1, 

and 

p{R,z) = (i?2 + £!)-7/2 ^ g-^r--'^(l-e%in^6l)-^/^ (4) 

The latter expression can be expanded in a power series in e^sin^^ using the binomial theorem, with the result: 

pir, e) = q<r-^ ^ 1 (7/2), (e^ sin^ e)\ (5) 

fe=0 

where {■ ■ ■)k is Pochhammer's symbol, which is defined in terms of Gamma-functions as {x)t = T{x + t)/T{x) (cf., e.g., 
Gradshteyn & Ryzhik 1994). 



2.2 Distribution functions 

2.2.1 Self-similarity 

DFs of the form f{£,L'^,L^) are sought, that generate the mass density in the potential (^. The integrals of motion are 
£ = — v'^ (the usual factor 1/2 in the kinetic energy term has been absorbed in the units), = r^{vg +v'^) and Lz = Rv^. 
We consider first the part of the DF that is even in Lz, fc{£,L^,L^). This part determines the mass density completely, 
because the latter is independent of a star's sense of rotation around the symmetry axis. 

The maximum angular momentum Lmax{£) at a given energy is attained by stars on circular orbits in the equatorial 
plane. The squared circular velocity is Vc{r) = /2, and hence 

r2 ... f i exp(-2f-l), 5 = 0; 

^.naxltj-| 1 [25£:/(2-5)](^"2)/5^ 0<5<1. 

Without loss of generality, the DF can be considered to be a function , i^^, j)^), with 

C' = LViLx(f), = L^JLl,J£), (7) 

so that < r?^ < < 1- 

The discussion can be restricted to DFs that are s elf- similar , since both the potential and the mass density are scale 
free. For such DFs there exist constants ci and C2 such that 



(8) 
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for all r, v and p. Following White (1985), we substitute fc{£ ,ri'^), differentiate with respect to p, and then set p = 1. This 
shows that ci = —5/2 for all < 5 < 1, and that /c must have the general form 

r /f ^2 2-, r^,^2 / exp(c2f ), 5 = 0; , . 

/.(£:,C ,r,) = F(C,,)x| ^^^^^^ (9) 

where F is an arbitrary non-negative function. 



2.2.2 DF components 

The mass density is the integral of the DF over velocity space, 



A^v f = I Avv'^ I dC sin^ / dr /o, (10) 

Jo Jo Jo 

where the variables {v,£,,t) are spherical coordinates in velocity space: 

Ur = wcos^, We = w sin ^ cos r, = w sin ^ sin r. (11) 

It is convenient to transform to the new integration variables 

£ = = 1 - (uV*), C^ssin^C, '7f = smT. (12) 
This results in 



p = *^/2/ A£{l-£y dCMl-C')"'^ / dr72(l-r72)-^''^(r72)-^''"/,. (13) 
Jipoo/* Jo Jo 

The lower limit of the outermost integral is ^'oo/^ = for < 5 < 1. For the logarithmic potential 5 = it is ^Poo/* = — oo. 
The integrals of motion {£,(^^,r]^) can be expressed in terms of the integration variables {£,(^^,1]'^) as follows: 

c ,T,T 2 . 2 a^^2 ^2 ,T, 272 /1 ■?^ . / 2 exp(2*£' + 1), 5 = 0; , , 

' /s , s s ^ ; I 2[25*£:/(2-5)]<2~*)/'', < 5 < 1. ^ ' 

We restrict ourselves in equation (^ to smooth functions F that can be expanded as sums of terms of the form C"^** rf^ . The 
entire DF is then a sum of self-similar components of the form 

/■[c2,M,A]/o ^2 2x^ >-2m 2A f exp(c2£:), 5 = 0; 

Substitution of equations (|l|) and (|l|) in equation (|l|), and carrying out the triple integration, shows that each component 
generates a mass density 



where the factors dI'^^.a'.a] given by 

£,[C2,M,A] ^ g( 1 ^ + 1 ) g( 1 1 _ ^ + ^) 



2^-" exp(A-M)r(|-/i-HA) (C2-Hf -2m+2A)-2+m-a^ ^ ^ Q; 

5-3/2 ^ 



-A 



(|)(2^)^ B(|-M+A,|[c2 + f-(2-5)(M-A)]-i), 0<5<1. 



(17) 



The function S(. ..,...) is the Beta-function, which is defined in terms of Gamma-functions as B{x, y) = T{x)T{y) /V{x -\- y) 
(cf., e.g., Gradshteyn & Ryzhik 1994). The 1)1=2. m, a] 

are continuous functions of 5 in the limit 5 [ Q. These results were 

obtained independently by Evans (1995, priv. comm.). 

Equation (|l6[) shows that in order to reproduce the mass density (H) with DFs of the form 



/e(£:,c',r?') = 5]^«'=^'''''i/^''''''(£,C',r,^), (18) 

one requires that 

C2=7-(35/2), and Y.a^^^'^''^^ D^^^'^''^^ ^ ^^{-f /2)u (for fc = 0, 1, 2, . . .). (19) 

In addition, the expansion coefficients qI'^^.m.a] jjj^g^ |-,g ^ero for all A 7^ fc (with fc = 0, 1, 2, . . .). The value of C2 ensures that 
the density components fall as r~^, for all < 5 < 1. Henceforth, C2 = 7— (35/2) is substituted in all equations that involve C2. 
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Table 1. Special cases of the function h{x^) defined in equation (26), and of the function j{x'^) for 5 = 1, defined in 
equation (32). The function [j(x^)]g—i for 7 = 4 can be reduced to an elementary expression only when 2/3 is an integer. 



4^3(1 2 , 2^ , 2^ ; 1 1 _ /3, 7 - /3 - i; x2) 



1 

1 I X arcsin x 

^ (l_^2).3/2 

1 + 



2-\-x^ I 3a: arcsin x 



1- 2/3+(l+2/3)3:^ 
(l-2/3)(l-a;2)2 

2- 2l3+(l+2fi)x^ 
2(l-/3)(l-x^)i! 

,2-2%lL2,-, -''-'£ {h£ [IS-'-'^Kx^) dx] } 



arcsin 3: 
2{l-,3){l-2:2)5/2 



2. 2. 3 Two families of DFs 

Many DFs can be constructed that satisfy equation (p^). Two particular cases are discussed here, which differ in the choice 
of components jJ^^.m,*;] rpj^^ fixsX set has ^ equal to the same constant for all components, so that /e is a separable function 
of E, L^/L^ax(-E') and L^/L^^{E). In the second set, the components are chosen such that fe is a separable function of E, 
/L^^^{E) and Lj./L^ . The latter models turn out to have a velocity distribution anisotropy that is independent of position 
(cf. Section 2.3 below). 



In case I the DF is built entirely with components /i' 
fl{£,e,n) = Cog{£)C^^ j{e^rf), 
where the functions j and g are defined as 



• / 2 2 



k=0 



/ 2 2\k 

ak{e Tj ) , 



exp{-y£), 

(5£)(7/*)-(3/2)^ 



for which ^ is equal to a constant p. The DF is then: 



5 = 0; 
0<5<1, 



and in addition Co = a''=2./3,o] ^^^^ ^ ^Ic2,i3.k] ^^^2k^ic2,H,o]y Comparison with equation shows that 
Co = qyD^'''^'°\ Ok = (7/2)fc ijI^^-'^'O'/ (fc! D^-2,fi,k]^^ 
where the factors d['^2,/3,'=] g^j-g given explicitly in equation ([17 



(20) 



(21) 



(22) 



In case II the DF is built entirely with components jJ'^^.M.fc] which fi — f3 + k, where /3 is a constant. The DF is then: 



where the constant Co and the function g[£) are as defined for case I, the function h is defined as 



(23) 



k=0 



bk (eW/e)", 



and bk =qI'=2./^+'=.'=)/(, 



j,[c2,/3 + fc,fe] 
£)[c2,/3,0] ^ 



,2k tc2, 13,0] 



). Equation ( |l9| ) now gives 



2Jk 



(24) 



(25) 



The second equality follows upon substitution of equation (^) and some algebraic manipulations. The series ( p^ ) is therefore 
recognized as a hypergeometric function: 



MeV/r) = 2i^i(l, 



2. 1.^2^2 



2 /a2 



Recall that the generalized hypergeometric function pFq is defined as 

(ai)fc • • • (Qp)fe x'' 



pFq (qi 



■ (Mk k\ ' 



(26) 



(27) 



for < x < 1, Q; > for / = 1, . . . ,p, Pm > for m = 1, . . . ,g. It sometimes reduces to an elementary function for special 
values of the parameters. One always has pFq (...;...; 0) = 1. The function h in equation ( p6| ) is elementary for integer values 
of 7. The explicit expressions for 7 = 1, 2, 3 and 4 are given in Table 1 (and are illustrated in Figure 1 below). 
For both case I and case II, the DF is positive definite if and only if both 



/3< 1, 



and 



7 > (5/2)+p{2-5). 



(28) 



If 7 > 3/2, then the latter constraint is satisfied for all /3 < 1 and < 5 < 1. The DF for either case is easily evaluated 
numerically using the series expansion of j or h, respectively. These series generally converge quickly. 
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2.2.4 Special cases 

There are special cases for which the case I and case II DFs simplify. Some are collected here. 

In the spherical case q = 1, one has j = h = 1, and fl and fl^ are identical. The DFs now depend only on £ and 
and not on L^. They describe constant-anisotropy models (e.g., Henon 1973). The value of (3 controls the velocity dispersion 
anisotropy (see Section 2.3). 

The (3 — > —00 limit yields the model with all stars on circular orbits (for which ("^ = 1): 



V 1- -fii tp\x /I , , 2 2^ _ V ^ / exp(7/2), 5 = 0; 

hm = hm = K-00 g{t) Su{l - C ) h(e rj ), K-00 = — x i (2-s\(s-i)/s ^ (29) 

/3^-oo TT^ (~2~j ' < < 1 

This model is physical for all relevant 7 and 5, cf. equation (^^. The function 5o{. . .) is Dirac's delta-function. The equality 
of the case I and case II DFs in this limit follows directly from the construction of the DFs. The case I DFs are built 
from components jJ'^^.A'.fc] ^j^jj ^ — ^j^g gg^gg jj ppg fj-om components with fi — P + k. For (3 — > —00 these components 
are identical. 

The /3 ^ 1 limit of the case II DF yields the model with all stars on radial orbits (for which (^^ =0): 

„7 f 2(7- 2)^/2 exp(l), (5 = 0; 

lim/," = 7^lg(£)&(C')(l-e^in^e)-^/^ Ki = ^ x l^,,,^^s ^^^j^^+s-^ 0<5<1 

^ t V2-4'' r(^[7+«-2]-i) ' - 

This model is physical only for 7 > 2 — 5/2, cf. equation (|2^). The quantity sin^ 9 is an integral of motion for radial orbits, 
and this DF is thus indeed a solution of the collisionless Boltzmann equation. The /3 — > 1 limit of the case I DF does not yield 
a model with all stars on radial orbits (see Figure 3 below). 

When /3 = 0, the case I DF depends only on £ and rj'^, and hence is independent of L^. This is the classical axisymmetric 
f{E,L,) model. 

When the potential is Keplerian [5 = 1), the function j that appears in the case I DFs can be expressed in terms of a 
generalized hypergeometric series. Equation (^) gives for this case 

^°-B(i,i)e(i,l-/3)S(|-;3,7-/3-i)' fc! (1)^(1-/3), (7-/3-1)^ ' ^'^^ ' ^ (31) 



The series ( |21| ) is thus 

j(eV)= 4^3(1, ^,2^f±i,^^f±2;i,l-/?, 7-/3- i;eV), (for 5 = 1). (32) 

This reduces to an elementary function when 7 and 2/3 are integers. Explicit expressions for 7 = 1,2,3 and 4 are given in 
Table 1 (and are illustrated in Figure 1 below). The expressions for 7 = 1, 2 and 3 are elementary for arbitrary (3. In the 
special case /3 = 0, one has (cf. eqs. |^ and ^]) that 

go= ^ 1 3. , i(eV) = 3F2(i,2±i,2±^;i,2a^;eV), (for 5 = 1, /3 = 0). (33) 

27rZ3(7-2, 2) 

This reproduces the f{E,Lz) scale-free large-radii limit given in equation (3.24) of Qian et al. (1995) for general 7, and in 
equation (B2) of Dehnen & Gerhard (1994) for 7 = 4. The elementary expressions for integer 7 follow from those given in 
Table 1 upon substitution of /3 = 0. 



2.2.5 Odd part 

The mass density determines uniquely the part of the DF even in L^. The odd part, fo{£,C,v), 

can be specified freely, with 

the only constraint that the total DF / = /c + /o be positive definite. A natural choice is 

U£,e, rj) = {2s - 1) sign(r?) r,^* U{£, C.rf), (34) 

where < s < 1 and t > Q are two free parameters. The fraction of stars on circular orbits in the equatorial plane that 
rotates clockwise is equal to s. A model with s = 1/2 is non-rotating. The parameter t determines the extent to which the 
net rotation of the model stems from high-angular momentum orbits. The odd part with f = and either s = or s = 1 is 
referred to as 'the maximally rotating odd part'. All stars with Lz 7^ have the same sense of rotation around the symmetry 
axis in a model with this odd part. The DF that generates the largest amount of mean streaming consistent with the given 
mass density and potential is the (3 ^ —00 model with the maximally rotating odd part. This model is referred to as 'the 
maximum streaming model'. 
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2.2.6 Density of states 

The amount of mass contributed by stars on orbits with given (£,("^,77^) is not determined solely by the DF, but also by the 
'density of states' (e.g., Binney & Tremaine 1987). In the present context the density of states w{£,(^^ ,1]^) is defined through 
the following expression for the total mass of the system: 

/■*{0) pi 

M= d£ AC dn'wi£,C,v')f{£,C\v^)- (35) 

J*^ Jo Jo 

To obtain an explicit expression for the density of states for the case of an axisymmetric mass density in a spherical 
potential, one expresses the total mass M as the integral of 27rr^ sm p{r, 6) over dr A9. Then p{r,0) is substituted from 
equation (p^), and the integration variables are transformed to {£■ ,C Rearrangement of the order of the integrations 
then yields an expression for the density of states as a two-dimensional integral over dr AO. For a spherical potential the 
integral over AO can be evaluated analytically, with the final result 

w{£,e,il^) = 2n' L^^4£) iCr'^^vT'^' [ ^'^'^ ' r {r^ [*(r) - £:]L-L(£:) - C'}"'^' dr. (36) 

Jr_(£,(^) 

The radii r± are the roots of the expression in curly braces. The interval [r_,r+] contains the radii accessible to a star with 
given {£,C). The integral in equation (^) can be calculated analytically only for the Kepler potential, 5 = 1. In this case 

w{£,e,v') - ^£''^'{CY'^'iv'r'^', (for 5 = 1). (37) 

Equation (^H]) can be used to calculate the differential mass distribution as function of the integrals of motion. For example, 
for a spherical mass density in a Kepler potential one has, for either the case I or the case II DF (cf. Section 2.2.3), / = 

^^^7-0/2) ^-2/3^ and thus 

^ = -5 = 1, g = 1). (38) 



2.3 Intrinsic velocity moments 

The intrinsic velocity moments {vl v]^ v"^) of arbitrary order follow from 

1 I m n\ /i3 J. I m n /nn\ 

p(VrVo v^) = Avfv^ve v^, (39) 

where 1,771,11 > are integers. The quantities p{u[ 17™ v^) with Z + m + n = 2 are often called the stresses. As before, we 
transform to the integration variables {£,(^'^,7]^) (eq. ^]), and use the relation 

vl Vl = $('+'"+")/2 (1 _ J)(!+™+n)/2 _ ^^^{m+n)/2 _ ^)™/2 (^)"/2, (40) 

which follows from equation (|ll|). For a DF component jl'^^,^,^] given in equation (p^, with £2=7— (35/2) as before, the 
integral yields 

p{vl v';)^"^'^'''^ = sin'^ 9 Z3(!2±i, A + ^) B(^, 



2 i^'i 2 -'^^2'-^ J^i^-i 2 

l + r,T.+ 



2^'^ exp(A-M) V{l-p + \+l^^^) (7-2At+2A)-2+"-^-^^^, 5 = 0; 



^-(3 + ;+m + n)/2 



2-5 ^ 



2 



A (41) 



so that the contributions of each DF component to the velocity moments are all elementary. 

The moments p{v\. v^) of the DFs with the case I and case II even parts defined in Section 2.2.3, and the odd part 
defined in equation (^^, can be expressed as series of these components: 

00 00 
p{vl v^f = (25 - 1) Co ^ e''= p{v\ w^)'^-"'^', p{v\ C = (25 - 1) Co ^ 6, 6^*= p{v\ ^;)[=-''+'='"i,(42) 

fc=0 fc=0 

where one should substitute: 5 = 1 and A = fc for even n; and S — s and \ — k + t for odd n. The summations are power 
series in sin^ 0. The velocity moments are easily evaluated numerically from these power series, which generally converge 
quickly. Substitution of the definitions of and bk in equation ( ji^ ) shows that the velocity moments of the case I and case II 
DFs are always identical on the symmetry axis. 
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Table 2. The intrinsic second order velocity moments. Here = sin^ d. The results are elementary for case II. For 
case I the second order velocity moments reduce to hypcrgcometric functions for the Kepler potential. These reduce to 
elementary functions when 7 and 2,8 are integers. This is illustrated by the results for 3 = listed in the bottom part 
of the Table. The function that occurs in the expression for {v^) for this case can be evaluated by means of the 

relation 3i^2(^, |; |, ^^■,x^) = £[x 2Fi(l, ^^■,x^)]. For 7 = 1 the listed expressions for {v^} and 

are independent of 13, and hence these results are valid for all /3. The quantity (i;^) does depend on /? for 7 = 1. Cases 

not covered in this Table can be calculated with the series expansion (42). 

Case II DF, arbitrary 5, P and 7 

{ (^-f ) ln(l x'^) (7 = 2) 

/ 2x n 2, / 2+ii^ln(l-x2) (7 = 2) 

{vD = (1 - X I ^ 



Case I DF, Kepler potential (5 = 1), arbitrary fS and 7 



p (2 7-2/3+1 . 7-2/3+3 2N 
\"r/ 2r(7-2/3+l) 2-^11,2' 2 ' 2 ' > 



2r(7 


-2/3+1) 


(1- 




2r(7 


-2/3+1) 


(1- 


x2)-r/2 



1 



(^i)= 2^F^TO (l-/3)4F3(l,t,^^,2-^;l-A2,2^;x^) 
("|>= felft^ (1-/3) 5F4(l,i,2^,|,2-^;i,l-A2,^^;x2) 

Case I DF, Kepler potential (5 = 1), f{E, L^) model (/3 = 0), and integer 7 

Iy2^ _ l2^ _ O^Z^Y!!! p,l 7±1.7±3. 2> / .2 x _ (^-f f 7 2±i 2.1 7±3. 2> 

\^T-/ — ~ 2t-(7+1) 2 1 2 / ~ 2r(7-i-l) 3-^2(,2i 2 ' 2' 2' 2 ' ^ / 

(l-a;2-|l/2 (^,^2., l_(l_^2)l/2 

2 1 ini+-) 1 (2-x2-i^lni±2) 

2rx-^ \ 2x 1 — x J 2rx^ \ x 1 — x J 

2 2-3x'^^x^-2(l-x^)^/^ -6+9a:^-2x^+6(l-x^)3/2 

4 (lz4)(3-2x2-3(l_:£!)lnl+2) (^(3x2 + -^-6+ 2(1^ In i+i) 

4t-x* \ 2a: 1 — x / 2rx^ \ l — x-^ x 1 — x } 

For the case II DFs the power series in sin^ B reduces as before to a generahzed hypergeometric function: 

X 3^2(1, f,r+^;i,l+r+21±21:e2 ^in^ 61) 

2 "xp(r 3 ^ — -^j— , (5 = 0; 

r(--»3) (^-2,3) C43'v 



X < 



j-(i + m+n)/2 



B(i-/3.7h-(2-5)/3]-i) 



< 5 < 1. 



For the case I DFs the power series in sin^ ^ reduces to a generahzed hypergeometric function only for the Kepler potential 
(5=1): 

p{vivTv;y = 

(2S l)r (4sm^)g ^^^^^^ e(|_^,^_^_i) 

X 7^6(1, 2, 2^f±i, 2^f±2,7-/3-i+T,T+^i±i, l-/3+r+2i±:i; i, 1-/3,7-/3-1, 

1 +r+ , T-f +^ +r+ i±ig±s , 7-2^/3+2 -I-T+ i±i|±s ; sin^ 6') , (for ,5 = 1). (44) 

In these equations one should substitute: 5=1 and T = for even n; and S = s and T = t for odd n. 

For the low-order velocity moments the pFq in these expressions often simplify. This is illustrated by Table 2, which lists 
second order velocity moments for the case I and case II DFs. For case II these are elementary for all 6, 7 and /3. For case I 
the second order velocity moments in a Kepler potential are elementary when 7 and 2/3 are integers. 

It is useful to consider Binney's (1980) velocity dispersion anisotropy function /3b, defined as 
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On the minor axis (sin 6* = 0) , and in the spherical hmit (e = 0) , one has 

{vl) = {vl), [3B{r,9)^p, (if sin' = 0), (46) 

for both case I and case II. For the special choice of case I and (3 = (i.e., the f{E, Lz) models) one has {vq) = {vl) everywhere, 
for both spherical and flattened models. More interestingly, it follows from the expressions in Table 2 that for case II, Binney's 
anisotropy function is independent of the polar angle: 

/3b(j-, 6*) = P (for case II and any e, 9). (47) 

This means that the case II DFs describe flattened constant-anisotropy models. Models with /3 = are isotropic, those with 
/3 < are tangentially anisotropic, and those with /3 > are radially anisotropic. 

2.4 Projected velocity moments 

Following Evans & de Zeeuw (1994), Cartesian coordinates {x',y',z') are defined such that x' lies along the projected major 
axis of the model, y' lies along the projected minor axis, and z' lies along the line of sight. The inclination of the galaxy is 
denoted by i, such that for i = 90° the object is edge-on, while for i = 0° it is face-on. 

The projected mass density Ep(a;',j/') on the plane of the sky can be calculated as in Fillmore (1986), which yields 

-1+1 

/oo / i2 \ — 2 — 

dz' p(T, 9) = ^Bm^,\)[x'^ . (48) 

-oo ?P V / 

The quantity Ep is proportional to the surface brightness, if the mass-to-light ratio is constant. The axial ratio of the similar 
concentric elliptical isophotes is q^, which satisfies = cos' i + sin' i. Their ellipticity is tp = 1 — gp. 
The line-of-sight velocity at any given point is v^i = C'rVr + Cgve + C^v^, where 

Cr = sin S cos sin i + cos 6' cos i, Ce = cos S cos (?!)sin i — sin cos i, C,^ = — sin (/!>sin i. (49) 

The n-th line-of-sight velocity moment {v"^,) at any given point satisfies 

p(<'> - E E ^ ^) ci cl c;~^'' piviv'ov;-'-'-), (50) 

j=o k=o V-'/ V / 

which is obtained by using the binomial theorem twice. The quantities p{v^V0V^~''~'') are as given in equation (^). The n-th 
projected line-of-sight velocity moment («")p on the {x' ,y') plane of the sky follows upon integration along the line of sight: 

{v")p = — / dz'p{?;",). (51) 

This integral must generally be evaluated numerically. 

2.5 Velocity profiles 

The velocity profile (VP) at any point {x',y') on the sky is the line-of-sight velocity distribution of the stars: 

VP{v,,) = ^ Jdz' J dv^, J dvy, f, (52) 

where / is the DF. The integration limits are set by the fact that the integrand / is zero for those points in phase space where 
f > The moments of the VP are equal to the projected line-of-sight velocity moments given in equation (^), i.e., 

dv,, [VP Mr = («?)p- (53) 

Observed VPs are often represented as a Gauss-Hermite series (e.g., van der Marel & Franx 1993; Gerhard 1993). This series 
is parametrized by the normalization 7g, mean V and dispersion a of the best- fitting Gaussian to the VP, and the Gauss- 
Hermite moments /i3, /14, . . ., which quantify deviations of the VP from this Gaussian. Calculation of these quantities for the 
models discussed here requires knowledge of the VP. The VP can in principle be calculated by direct numerical evaluation of 
the triple integral (^^, but this is cpu-intensive and not convenient in practice. An alternative is to recover the VP from its 
moments {v"i)p, which can be obtained by single integrations (cf. eq. [Q). This is a well-known, ill-conditioned mathematical 
problem, but after some experimenting an algorithm was found that is satisfactory for the purpose at hand. It is described 
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logio j for |S = -3 logio j for |S = 




Figure 1. The functions j(e^rf) and h{e'^ rj^ / C,'^) responsible for the flattening of the mass density for the case I and case II DFs, 
respectively. Solid curves are for the Kepler potential (5 = 1), dashed curves for the logarithmic potential (<5 = 0). Each panel has four 
curves for each choice of potential, corresponding to mass density power-law slopes 7 = 1, 2, 3 and 4. The function j depends on /3. It 
is shown for three representative values: /3 = —3, and 0.75. The 7 = 1 curves are absent in the /3 = 0.75 panel, because these do not 
correspond to a physical DF. The function h is independent of both /3 and 5. The vertical scale in all panels is logarithmic. 

in Appendix A. It works well, except for the small region of parameter space describing strongly flattened models with a 
logarithmic potential and large anisotropy, which will not be discussed in the remainder of the paper. 



3 MODEL PROPERTIES 
3.1 Distribution functions 

In Section 2 two families of DFs were presented (referred to as case I and case II) which generate a scale-free spheroidal mass 
density with power-law slope 7 > and flattening 5 < 1, in a scale- free spherical potential with power-law slope < 5 < 1. 
The part of the DF even in has (for each family) one free parameter ~oo < /3 < 1, which regulates the dynamical structure 
of the model. The DFs are given by equations ( po| ) and (p^), respectively. They are physical if and only if equation (^^ is 
satisfied. A convenient ad hoc choice for the odd part of the DF is given by equation (p^. This odd part has additional free 
parameters s and t, which regulate the mean azimuthal streaming in the model. 
The DFs fl and have a factor CqC'^'^ 9{£) 

in common. The normalization constant Co depends on 7, 5, P and q. 
The quantity (^^ is defined as the ratio L'^ / L'^^^{£) (cf. eq. Q). The function g{£) is a scale-free function of the energy £, as 
required by the nature of the density and potential. It is fully determined by 7 and 5. 

In the spherical case one has fl = /e" = CoC"^'^5(5)- For flattened models the case I DFs have an extra factor j{e^rf), 
while the case II DFs have an extra factor hie^rf /C,"^). The quantity rf is defined as Ll/L'^^^{£) (cf. eq. |Q]). Note that 
jC^ _ 7;,2^/^2 rpj^g functions j and /i, respectively, are responsible for the flattening of the mass density. The axial ratio 
enters into these functions only through the eccentricity e in the argument. The function j (case I) depends on 7, /3 and 5. 
The function h (case II) depends only on 7. 

Figure 1 displays the functions j and h for 5 = and 1, and 7 = 1, 2, 3 and 4. The function j is shown for three 
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Figure 2. Velocity moments for a model with a Kepler potential (<5 = 1) and a mass density with flattening q = 0.8 and power-law 
slope 7 = 4. Shown are: (i)^)-*-/^ (solid curves), (vg)^^^ (dotted curves), {v'^)^^^ (short-dashed curves) and (ti0)max (long-dashed curves) 
as function of polar angle 8 at radius r = 1. Top panels are for the case I DFs, bottom panels are for the case II DFs. From left to right 
the value of the model parameter /3 is: —3, and 0.75, respectively. In the top middle panel, when / = f{E, Lz), one has (v^) = {vg). 
The symmetry axis has 6 = 0°, the equatorial plane has 6 = 90°. 



representative values of /3: —3, and 0.75. Both j and h increase monotonically as function of their argument. At fixed 
flattening, the physical range of the argument runs from to e'^ = 1 — g^. For realistic elliptical galaxy models {q ^ 0.3), the 
functions j and h can vary by as much as two to three orders of magnitude over their physical range. However, this does not 
imply that most of the stars in the system are on orbits with either rj^ = 1 (case I) or rj^ = (^^ (case II), respectively, because 
the density of states for these orbits is low (cf. eq. |36| ). 



3.2 Intrinsic velocity moments 

To understand the dynamical structure of the models it is useful to focus on the first and second intrinsic velocity moments. 
These are easily calculated for any combination of model parameters using the formulae in Section 2.3. As an example 
consider the particular case 7 = 4, g = 0.8 and 5 = 1. Figure 2 shows for the case I and case II DFs, for P = —3, and 
0.75, the dependence of (v^)^^^, {vg}^''^ and (v^)^''^ on the polar angle 6, at radius r = 1. The dependence on r is simple 
(see Table 2), because of the scale-free nature of the models. The mean azimuthal velocity («^)max for the maximally rotating 
model associated with this even part is also shown in the figure. 

The case II DFs have a constant ratio of rms radial to rms tangential {v^ = vg + v'^) motion as function of 6. This ratio 
is determined by the model parameter (3, cf. equation (|47|). The models with (5^1 have only radial orbits (with = 0), 
while the models with (3 — > —00 have only tangential orbits (with = L^^^{£)). The quantity (v^) is constant as function 
of 9 for the case II DFs, and hence so is {vg + v^). What does vary as function of 6 is the ratio {v^) /{vg). It is unity on the 
minor axis, and increases monotonically with 9. 

On the symmetry axis the velocity moments for the case I DFs are identical to those for the case II DFs. Away from 
the symmetry axis they behave differently. The ratio {v^)/{vg) for case I increases monotonically with 9, as for case II. By 
contrast, the ratio {v^) / {v^) also increases monotonically with 9, rather than being constant, as for case II. This is not very 
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case I case II 




/3/(2-/3) iS/(2-/5) 



Figure 3. Ratios of various velocity moments as function of /3/(2 — /3), for a model with a Kepler potential (5 = 1) and a mass density 
with flattening q = 0.8 and power-law slope 7 = 4. The left edge of each panel corresponds to /3 ^ -co, the right edge to /3 — > 1. Shown 
are the values in the equatorial plane of [(?-'^)/ (^r)]^'^'^ (solid curves), [(^^)/ (^^)]^^^ (dotted curves), (^,^)max /{Vr+v^y^^ (short-dashed 
curves), and also the value of = [{'"g) / {'"r)]^^'^ the symmetry axis (long-dashed curves). The case I DFs are always 

more tangentially anisotropic in the equatorial plane than the case II DFs, while on the symmetry axis they are identical. 



pronounced in the top left panel of Figure 2 for /3 = —3, because the case I and case II DFs are identical in the limit /3 — > — oo. 
However, it is very clear in the top right panel for /3 = 0.75. In fact, the case I /3 = 0.75 model has (Vf) < (v^) on the symmetry 
axis, and (vt) > {v^) in the equatorial plane. 

The dependence of the intrinsic first and second order velocity moments on /3 is further illustrated in Figure 3. This figure 
shows ratios of various velocity moments in the equatorial plane and on the symmetry axis, as function of /3/(2 — j3) (this 
choice of abscissa is useful because it maps the infinite interval — oo < /3 < 1 to the finite interval —1 < /3/(2 — /3) < 1). The 
figure clearly demonstrates the equality of the case I and case II DFs for (3 — oo. It also shows that the case I DFs are always 
more tangentially anisotropic in the equatorial plane than the case II DFs, while on the symmetry axis they are identical. 
The ratio {u</))max/{'u? + v'^)^^'^ of the mean azimuthal streaming and the rms motion in the equatorial plane {6 = 7r/2) is a 
monotonically decreasing function of jS. The maximum possible relative importance of mean streaming thus decreases as the 
importance of radial pressure in supporting the shape of the system increases. 



3.3 Projected velocity moments 

A useful observational indicator of the dynamical structure of a stellar system is the ratio v of the rms projected line-of-sight 
velocity on the major and minor axes (van der Marel 1991): 

U = (u^/)p,maj / (-uf')?,™!!! = (o^p + Up)maj / (o-p)min, (54) 

where {v1i)-p is defined in equation (|5^), and Vp and CTp are the observed mean streaming and dispersion. This ratio depends 
only on the even part of the DF. It is generally a function of radius. However, in our scale-free models it has a constant value, 
which can be evaluated numerically as described in Section 2.4. 

Figure 4 shows v as function of /3/(2 — 13) for the case I and case II DFs, for a system with 7 = 4, (5 = 1 and projected 
axial ratio = 0.8. The different curves correspond to different values of the intrinsic axial ratio q, and hence to different 
inclination angles. For the case I DFs v is generally an increasing function of /3, although v is close to constant for —00 < /3 0. 
The flatter case I models with smaller inclination angles have larger v. For the case II DFs is a decreasing function of /3, 
the more steeply so for the flatter models with smaller inclination angles. The results in Figure 4 are generic for other values 
of 7, 5 and q^. 

Following Binney, Davies & lUingworth (1990), van der Marel (1991) used solutions of the Jeans equations to compare 
the predictions of models with / = f{E,Lz) to kinematical data for 37 elliptical galaxies. He concluded that these models 
generally predict values of 1/ that are too large compared to the observed values, the more so for smaller inclinations. For 
the models discussed here, / = f{E,Lz) corresponds to case I with /3 = 0. Figure 4 shows that none of the case I models 
can produce values of v that are appreciably smaller than those of the edge-on f{E,Lz) models. From this it follows that 
the case I DFs (or their self-consistent generalizations) are probably not a good representation of real elliptical galaxies. This 
can be attributed to their property that the ratio of rms tangential to rms radial motion always increases strongly with 6. 
Apparently, this is not realized in nature, although it does lead to dynamically acceptable models. The galaxies in the van 
der Marel (1991) sample have values of v roughly between 0.9 and 1.3. The case II DFs in Figure 4 can easily reproduce this 
range. If elliptical galaxies have DFs similar to that of the case II models, Figure 4 indicates that they will most likely have 
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Figure 4. The ratio u of the rms projected Hne-of-sight velocity on the major and minor axes, as function of /3/(2 — for a model with 
a Kepler potential (5 = 1) and a mass density with power-law slope 7 = 4 and projected axial ratio qp = 0.8. The left edge of each panel 
corresponds to /3 — > — oo, the right edge to /3 — > 1. The curves correspond to values of the intrinsic axial ratio q = 0.3, 0.4, 0.5, 0.6, 0.7 
and 0.8, each corresponding to a model viewed at a different inclination angle. The dotted lines correspond to u = I. In real galaxies u 
is generally between 0.9 and 1.3, indicating they are probably best fit by case II DFs with /3 > 0. 



case I, qp=0.8, max. rotation case II, qp=0.8, max. rotation 




P/{2-l3) P/{2-l3) 



Figure 5. Ratio Vp/up of the mean streaming and velocity dispersion on the projected major axis as function of /3/(2 — for the same 
set of models as in Figure 4, using the maximally rotating odd part. The left edge of each panel corresponds to /3 — * — oo, the right edge 
to /3 — > 1. Less streaming is possible in models with large /3, i.e., in models with more radial motion. 



/3 <: 0. This is consistent with expectation based on N-body simulations of galaxy formation through dissipationless collapse 
(van Albada 1982; Bertin & Stiavelh 1993). 

Figure 5 shows the ratio Vp/ap of the mean streaming and velocity dispersion on the projected major axis for the same 
set of models as in Figure 4, for the maximally rotating odd part. Models with lower inclination and smaller axial ratio q 
generally have larger Vp/ap, in spite of the fact that they have less of their intrinsic streaming along the line of sight. For 
both DF families Vp/ap decreases with /3. In models with more radial motion one thus expects to see relatively less streaming. 
Bright elliptical galaxies generally have Vp/ap ^ 0.4. Figure 5 therefore shows that bright elliptical galaxies rotate much slower 
than allowed dynamically, as is well-known (e.g., Binney 1976). 



3.4 Velocity profiles 

The Gauss-Hermite coefficients that characterize the VP shapes of the models can be calculated as described in Section 2.5 
and Appendix A. As an example, consider the model with 7 = 4, 5 = 1, axial ratio q = 0.8 and inclination i = 90°. Based on 
the results of the previous section, the discussion is restricted to the case II DFs. The parameter s of the odd part is varied 
from 5 to 1, while t is set to zero (i.e., /o is equal to times a step function, cf. eq. jsij). This yields models that range 
from non-rotating to maximally rotating. Figure 6 shows the Gauss-Hermite coefficients hg and h4 for different values of 
p. The abscissa in the figure is the observationally accessible quantity Vp/ap, which increases monotonically with the model 
parameter s. The predicted Gauss-Hermite coefficients depend in a complicated way on the model parameters 7, S, q and i, 
but none the less, the results in Figure 6 are generic for a wide variety of parameter combinations. 
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Figure 6. The Gauss— Hermite coefRcients and /14 on the projected major axis, as function of the ratio Vp/ap of the mean streaming 
and velocity dispersion. These coefficients measure deviations of the VPs from a Gaussian. The model has a Kepler potential (5 = 1), a 
mass density with power-law slope 7 = 4 and axial ratio q = 0.8, and is viewed edge— on. The curves arc for models with a case II DF, 
with /3 = -3.00, -1.83, -1.00, -0.41, (dashed curves), 0.29, 0.50, 0.65 and 0.75, respectively. The odd part of the DF has t = 0, while 
the parameter s is varied to produce models with different Vp/ap, ranging from non-rotating to maximally rotating. 

As argued in Section 3.3, the models that best fit real galaxies are probably those with /3 ^ in which the rotation is 
significantly less than the maximum possible. Figure 6 shows that these models predict —0.1 ^ h4 ^ 0.1. Opposite signs are 
predicted for hs and Vp/op, provided that /3 is not too close to unity. These predictions agree well with the observations of 
nearly all galaxies for which VP information is available (e.g., van der Marel & Franx 1993; van der Marel et al. 1994a; Bender 
at al. 1994). 

The even part of the VP is fully determined by the even part of the DF, and hence is independent of either s, t or Vp/up. 
The fourth-order Gauss-Hermite moment of this even part is sometimes referred to as Zi (e.g., van der Marel et al. 1994b). 
In Figure 6 its value is read off as the value of 7i4 at Vp/ap = 0. 



4 AN APPLICATION 

As an application of the models, consider the issue of the dynamical detection of dark matter in elliptical galaxies. Both 
tangential anisotropy and the presence of a dark halo can cause the observed velocity dispersion (jp to remain roughly 
constant as function of galactocentric distance R' , out to well beyond the effective radius R'^g. Proving the presence of a dark 
halo therefore requires the construction of anisotropic axisymmetric models, to rule out the possibility of strong tangential 
anisotropy. 

Carollo et al. (1995) presented stellar kinematical data for the four elliptical galax;ies NGC 2434, 2663, 3706 and 5018, 
out to ~ 2J?gff . Here the discussion is restricted to two of these galaxies, NGC 2434 and NGC 3706. Carollo et al. interpreted 
their data by constructing flattened models with / — f{E,Lz). From combined modelling of the major axis rms projected 
line-of-sight velocity (vl,)p = ap + Vp, and the Gauss-Hermite coefficient Z4 they concluded that the data for neither galaxy 
can be fit by any DF without invoking the presence of a massive dark halo. They also showed that the dark halos must be 
flattened, if the observed VPs are to be fit with an / = f{E, L^) DF. 

The families of DFs presented here can be used to further interpret the observations for NGC 2434 and NGC 3706. 
The model parameter 7 is chosen to fit the observed surface brightness slope at large radii, and qp is chosen equal to the 
average apparent flattening of the isophotes outside half the effective radius. This yields 7 = 2.94 and qp — 0.92 for NGC 
2434, and 7 = 3.36 and qp = 0.65 for NGC 3706. The potential is chosen to be logarithmic {5 = 0), since it has already 
been demonstrated that both galaxies must be embedded in a dark halo. This yields a flat ap proflle. Following Carollo et 
al. we study the Gauss-Hermite coefficient 2:4. This quantity is fully determined by the even part of the DF, and the only free 
parameters of the models are thus /3 and the inclination angle i. 

Figure 7 shows the model predictions as function of the assumed inclination. The data points at large radii fall between 
the two horizontal dotted lines (these represent the error bars at the outermost measured points, and are a conservative 
estimate of the observational uncertainty in the mean Z4 for R' > R'i.g; additional systematic errors in the spectral analysis 
due to template mismatching or continuum subtraction are believed to be at most |A24| ^ 0.03). The dashed curve shows 
the predictions for the f{E,Lz) model, i.e., case 1 and l3 — 0. The predicted 24 values fall below the observations for both 
galaxies, as shown already by Carollo et al. They demonstrated in addition that f{E, L^) models with a flattened dark halo 
do predict the correct za for both galaxies. The required flattening of the dark matter density is g « 0.7 for NGC 2434 and, 
somewhat implausibly, q ^ 0.3 for NGC 3706. The results of Figure 7 show that the data can also be fit with a spherical 
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Figure 7. The Gauss-Hermite coefficient Z4, of tlie major axis VPs of NGC 2434 and NGC 3706. This coefficient quantifies the lowest- 
order deviations of the even part of the VP from a Gaussian. A value Z4 < indicates that this even part is more flat-topped than a 
Gaussian, and a value 24 > that it is more centrally peaked. The data points outside the eflfective radius fall between the horizontal 
dotted lines. The curves show the model predictions as function of the assumed inclination. The lowest inclination plotted corresponds 

to an intrinsic axial ratio q = 0.3. Dashed curves are for f{E, Lz) models, i.e., case I and /3 = 0. These models do not fit the data. Solid 
curves are for case II DFs and various values of /3 (/3 = —3.00, —1.83, —1.00, —0.41, and 0.29, and in the right panel also (3 = 0.5). For 
these models a range of inclinations and (3 values provides an acceptable fit to the data. Our algorithm for calculating 24 (Appendix A) 
does not work well for very large values of /3 in a logarithmic potential, but calculations in other potentials indicate that Z4, keeps 
increasing monotonically with fi for ^ — ► 1. 



dark halo and a more general DF. The solid curves show the predictions for the case II DFs with different values of (3. For 
both galaxies there is a significant range of inclinations and values of /? for which the observed Zi values are matched. Hence, 
the VP shapes can be fit with both flattened and spherical dark halos. More data, e.g., along the minor axis, are required to 
distinguish between the various possible models. 



5 DISCUSSION AND CONCLUSIONS 

A study has been presented of stellar dynamical models for scale-free flattened spheroidal mass densities in scale-free spherical 
potentials. The mass density is characterized by its power-law slope 7 and its flattening q, while the potential is characterized 
by its power-law slope 6. The general form of the DFs was derived, and two particular families of DFs were studied in 
detail. The DFs of these families are separable functions of the integrals of motion, or combinations thereof. Both families 
have a free parameter — cxd < 13 < 1 which regulates the velocity dispersion anisotropy. In the spherical limit they reduce 
to constant-anisotropy models of the type discussed by Henon (1973). The DFs of the models can be expressed in terms of 
generalized hypergeometric functions or power series with known coefficients, which reduce to elementary functions in many 
cases of interest (Tables 1 and 2). Because of their simple structure, it is relatively straightforward to calculate the intrinsic 
and projected velocity moments. The latter can be used to reconstruct the projected VP shapes. 

The 'case I' distribution functions are anisotropic generalizations of the flattened f{E, Lz) model, which they include as 
a special case. The 'case IF distribution functions generate flattened constant-anisotropy models. For both families, Binney's 
function /3b on the minor axis is equal to the DF parameter (3. For the case I DFs, the ratio {vt)/{vr) increases monotonically 
with 0. For the case II DFs it is independent of 0, and the function /3b is everywhere equal to the parameter (3. The case I 
and case II DFs are identical in the limit (3 —00, which is the model built exclusively with circular orbits. Case I DFs with 
(3 = correspond to the classical f{E,Lz) DF. Case II DFs with /3 — + 1 correspond to models with all stars on radial orbits. 
Calculations of the ratio of the rms projected velocity on the projected major and minor axes show that real elliptical galaxies 
are probably best described by the case II DFs with /3 ^ 0. Such models also predict VP shapes consistent with observations. 

Two important conclusions can be drawn: (i) flattened axisymmetric stellar systems can have a large range of physical 
DFs and dynamical structures; and (ii) only a small subset of the possible dynamical structures is realized in nature. This 
agrees with the work of Dehnen & Gerhard (1993), who constructed self-consistent three- integral DFs for a flattened isochrone 
model. They restricted themselves to 'quasi-separable' functions of the integrals of motion, while the present paper has been 
restricted to two special families of DFs. The full range of possible DFs for flattened systems is therefore even larger than 
that discussed in either paper. 

As an application, the models are used to interpret the VP data obtained recently by CaroUo et al. (1995). They showed 
that the galaxies NGC 2434 and NGC 3706 must have dark halos, and that the dark halos must be flattened if the observed 
VPs are to be flt with an / = f(E,Lz) DF. Our models demonstrate that the data can be fit equally well with a spherical 
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dark halo, provided that the DF is more general than / = f{E, Lz). In particular, the case 11 DFs with (3^0 provide a good 
fit. Data along more position angles are required to discriminate between the various possible models and dark halo shapes. 

A disadvantage of the models discussed here is that the potentials of real galaxies are generally not spherical, especially 
not in the inner regions, where the potential is dominated by the luminous matter. This introduces a number of systematic 
differences with respect to the predictions of self-consistent models. For example, the velocity ellipsoids of the models always 
align with spherical coordinates, whereas this need not be the case in self-consistent fiattened models, although it often is 
a good approximation (Dehnen & Gerhard 1993; de Zeeuw, Evans & Schwarzschild 1995). Also, the tensor virial theorem 
dictates that a flattened mass density in a flattened potential has more rms motion parallel to the equatorial plane than the 
same mass density in a spherical potential. On the other hand, the potentials generated by flattened mass distributions are 
always more nearly spherical than the density distribution itself, especially at large radii where the monopole component of 
the potential dominates. Indeed, the models illustrated in Figures 1 to 6 (Kepler potential, mass density power-law slope 
7 = 4 and flattening q — 0.8) are the asymptotic large-radii limit of the self-consistent models studied by Dehnen & Gerhard 
(1993). The potentials of real galaxies are probably dominated by dark halos, at least in the outer parts, and these may well 
be nearly spherical. 

One case where the present models are certainly applicable is to the central density cusp structure around a nuclear black 
hole, where the potential is known to be spherical and Keplerian. From equation ( psj ) with (5 = 1 it follows that the physical 
DFs must have /3 < 7 — ^, for either the case I or case II. So at least for the particular families of DFs studied here, the 
presence of a shallow density cusp (7 < |) around a central black hole, precludes a large radial velocity dispersion anisotropy. 
Conversely, the constraint on /3 implies that oblate f{E, L^) models (i.e., case I, /3 = 0) are physical only if 7 > |, as shown 
previously by Qian et al. (1995). 

A useful generalization of the present work would be to build triaxial models in spherical potentials. This can be achieved 
by using DF components that involve powers of L^, Ly and L1, rather than just and (Mathieu et al. 1995; Evans, 
private communication) . 
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APPENDIX A: RECONSTRUCTION OF A VELOCITY PROFILE FROM ITS MOMENTS 

Various approaches exist for the recovery of a VP (^) from its moments (^l|). One possibility is to reconstruct the VP from 
its Gram-Charlier series expansion (e.g., van der Marel & Franx 1993; Magorrian & Binney 1994). This works well, but only 
for VPs that are close to a Gaussian. In the present context it did not always produce satisfactory results. An alternative 
approach was therefore used, in which the VP is represented by its value on a discrete array of N velocities. The VP values 
are then estimated from the first N velocity moments upon inversion of the so-called VanderMonde matrix (e.g., Press et 
al. 1992, section 2.8). The observable quantities are easily calculated from the discretized VP. This process was performed 
iteratively, using higher and higher A'^, until convergence was achieved in (70, V, cr, /13, /14). This generally works well, despite 
the ill-conditioned nature of VanderMonde matrix inversion. The reason for this is that the Gauss-Hermite series bears a 
strong resemblance to a Fourier series. Gauss-Hermite moments hi of higher order measure power in the VP on higher 
frequencies (Gerhard 1993). For high N, numerical errors cause the VanderMonde solution to become oscillatory on the scale 
of the velocity array spacing. This, however, does not influence the observables {•ycV, a, h^, h4), which only measure power 
on relatively low frequencies. 

The VanderMonde algorithm' was tested in various ways. A program was written that calculates the VP for the spherical 
(g = 1) case of the models by direct evaluation of the three-dimensional integral in equation (js^). The (70, V, cr, /13, ^14) 
calculated with the VanderMonde algorithm were found to be accurate, except for the case of a logarithmic potential and 
large anisotropy (/3 ^ —4.0 or /3 ^ 0.8). This has two reasons: (i) the logarithmic potential has no escape velocity so that it 
is more difficult to represent the VP on a finite velocity array; and (ii) the VP becomes discontinuous in the limit /3 ^ —00 
(only circular orbits), and singular in the limit (3 = 1 (only radial orbits) (e.g., van der Marel & Franx 1993). A second test 
was provided by the f{E,Lz) case (/3 = 0), for which the three-dimensional VP integral (|5^) could be calculated as in Qian 
et al. (1995). Again, the {'yG,V,a,hz,hi) calculated with the VanderMonde algorithm were generally found to be accurate, 
with the exception of very flattened models q ^ 0.6 in a logarithmic potential. This is understood from the fact that the 
VPs of these models become contrived and double-peaked for large flattening (Dehnen & Gerhard 1994). In both test cases 
inaccurate results were accompanied by non-zero values of h\ and /12, which should be zero by definition. For the general case 
we therefore took the values of hi and /12 as an indicator of the numerical accuracy of the VanderMonde algorithm. From 
this it was found that the algorithm fails only for the case of a logarithmic potential, large flattening and large anisotropy. 
No attempt was made to develop an algorithm to recover the VP in a more sophisticated way. This probably requires some 
form of regularization of the problem (see, e.g., Merritt & Tremblay 1994), which is outside the scope of this paper. 

This paper has been produced using the Royal Astronomical Society/Blackwell Science style file. 



